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1  SUMMARY 


The  goal  of  this  project  is  the  development  of  sophisticated  road  map  extraction  and  fusion 
algorithms  which  use  the  ground  target  radar  tracks  as  building  blocks.  We  develop  inno¬ 
vative  frameworks,  methodologies  and  algorithms  that  intelligently  fuse  tracks  of  different 
targets  which  move  along  the  same  road  over  time,  to  obtain  accurate  and  refined  estimate 
of  the  road  map. 

In  this  project,  modeling  the  roads  as  piecewise  linear  segments,  an  appropriate  data 
structure  is  developed  to  represent  the  road  network  map.  An  maximum  likelihood  (ML) 
road  estimator  using  radar  tracks,  as  well  as  its  corresponding  Cramer-Rao  lower  bound 
(CRLB),  is  derived.  To  properly  fuse  target  tracks  with  the  road  network,  a  track  to  road 
correlation/association  approach  is  proposed  to  align  track  segments  to  existing  road  seg¬ 
ments.  To  enhance  the  association  performance,  statistics  and  information  extracted  from 
the  track  estimate,  including  both  position  and  velocity  estimates,  are  used  in  the  track  to 
road  association  approach.  If  the  track  has  overlap  with  any  existing  road  segments,  a  track 
to  road  fusion  algorithm  is  proposed  to  fuse  the  track  with  the  associated  road  segments 
to  improve  the  road  map  accuracy;  otherwise,  a  new  road  is  added  to  the  road  network 
based  on  the  target  track.  This  process  is  repeated  till  all  the  available  radar  tracks  are 
integrated/fused  with  the  radar  map.  In  this  project,  we  derive  theoretical  results  for  road 
estimation  based  on  radar  tracks,  track  to  road  association,  road  to  road  association,  and 
road  fusion.  Numerical  results  are  provided  to  demonstrate  the  effectiveness  of  the  proposed 
road  estimation  and  fusion  approaches. 

2  INTRODUCTION 

In  this  project,  we  seek  to  extract  road  maps  based  on  ground  vehicle  tracks  which  are  es¬ 
timated  using  radar  measurements,  such  as  those  obtained  from  the  ground  moving  target 
indicator  (GMTI).  Accurate  and  up-to-date  roadmaps  are  crucial  for  many  purposes  such 
as  navigation,  target  tracking  [1,2],  and  airborne  knowledge-aided  space-time  adaptive  pro¬ 
cessing  (STAP)  [3,4],  Digital  road  maps  produced  by  the  National  Imagery  and  Mapping 
Agency  (NIMA)  as  well  as  the  United  States  Geological  Survey  (USGS),  the  two  most  com¬ 
mon  sources  of  such  maps,  often  have  errors  that  are  large  compared  to  the  resolution  of 
the  GMTI  sensors  [5].  Road  maps  may  be  manually  extracted  from  aerial  photographs  or 
synthetic  aperture  radar  (SAR)  images,  but  this  process  is  extremely  slow  and  laborious. 
Automatic  transformation  of  images  into  digital  road  maps  is  possible,  but  is  a  difficult  task 
and  prone  to  errors,  such  as  the  false  roads  caused  by  boundaries  in  the  image.  As  a  result, 
radar  track-driven  road  map  extraction  is  an  important  and  challenging  problem  in  scenarios 
where  road  maps  are  not  available,  or  where  they  are  not  sufficiently  accurate  or  up-to-date. 
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2.1  Literature  Review 

2.1.1  Road  Map  Extraction  from  GMTI  Data 

A  literature  search  has  been  conducted  for  road  map  conflation  using  radar  data.  There  are 
a  few  publications  on  the  topic  of  road  map  extraction  based  on  GMTI  data.  In  [5],  O’Neil 
modeled  the  road  itself  as  a  trajectory  through  space,  which  is  indexed  by  arc-length  instead 
of  by  time.  An  iterative  procedure  to  estimate  the  road  using  GMTI  data  was  proposed, 
which  starts  with  the  initial  knowledge  of  the  road  trajectory  estimate  provided  by  target 
tracking  algorithms. 

Koch  et  al.  [6]  applied  track  smoothing  techniques  to  improve  the  target  trajectory  es¬ 
timates,  which  are  used  to  build  the  road  map.  The  road  map  information  was  further 
employed  to  assist  tracking  of  targets  that  move  along  the  roads. 

In  [7],  Sklarz  et  al.  deemed  track  and  the  road  as  composite  curve  entities  rather  than 
sequences  of  points,  which  were  associated  and  fused  as  a  whole.  In  [7],  by  using  dynamical 
time  warping  (DTW),  similarity  measures  are  proposed  to  measure  the  distance  between  two 
curves,  which  facilitate  the  track  to  road  association  process.  The  track  to  road  fusion  prob¬ 
lem  was  cast  as  a  curve  fusion  problem  and  a  preliminary  curve  to  curve  fusion  approach  was 
developed  by  simply  fusing  matched  points  and  ignoring  the  unknown  correlation  between 
the  road  estimation  error  and  the  track  estimation  error. 

2.1.2  Road  Map  Extraction  from  Images 

Some  recently  published  papers  are  related  to  the  road  network  conflation  problem  using 
image  data  [8,9].  In  [8],  the  authors  proposed  to  extract  road  intersections  and  terminations 
from  imagery  based  on  spatial  contextual  measures.  The  extracted  intersections  were  used 
as  control  points  for  conflation  algorithm.  In  [9],  an  EKF  has  been  combined  with  a  special 
particle  filter  to  regain  the  trace  of  the  road  beyond  obstacles,  as  well  as  to  find  and  follow 
different  road  branches  after  reaching  to  a  road  junction. 

2.2  Background  in  Track  to  Road  Association 

Track  to  road  association  is  a  very  crucial  step  in  the  road  map  generation.  Here  we  first  study 
the  hypothesis  testing  problem  to  determine  whether  the  target  state  estimate  at  a  particular 
time  belongs  to  an  exiting  road  segment  or  not.  Below,  we  present  three  approaches  based 
on  the  normalized  distance  (Mahalanobis  distance),  likelihood  with  diffuse  prior,  and  the 
generalized  likelihood  and  discuss  the  relationship  between  the  first  two. 

2.2.1  Traditional  Track-to- Track  Association 

A  straightforward  way  for  traek-to-road  association  (correlation)  is  to  utilize  the  position 
estimates  contained  in  the  target  tracks.  More  specifically,  similar  to  the  track  to  track 
association  approach  in  [10],  we  propose  the  following  hypothesis  test  to  decide  whether 


Approved  for  Public  Release;  Distribution  Unlimited. 
2 


a  target  at  a  certain  time  is  on  the  existing  roads  or  not.  First,  let  x&  denote  the  target 
position  estimate  at  time  k  with  a  covariance  matrix  Pfc,  and  x;  denote  the  position  estimate 
of  an  arbitrary  point  on  the  existing  roads  in  the  road  map  with  a  covariance  matrix  P/. 
For  simplicity,  we  assume  that  x*.  and  x;  are  jointly  Gaussian  with  cross-covariance  matrices 
P =  P[k.  We  define  the  difference  between  these  two  position  estimates  as 

A  ki  =  xfc  -  X;  (1) 

If  both  Xfc  and  x;  are  unbiased  estimates  and  they  correspond  to  the  same  location,  then 
A M  follows  a  Gaussian  distribution  with  zero  mean  and  covariance  matrix  Pfc  +  P;  —  P^- P;*.. 
The  point  in  the  target  track  at  time  k  can  be  associated  to  the  existing  roads  by  testing 
the  following  statistic: 

di  =  min  [Ay  (Pfc  +  P;  -  P h  -  P ik)~l  AfcZ]  (2) 

Note  that  the  search  space  of  the  above  optimization  problem  can  be  significantly  reduced,  if 
a  gating  procedure  is  used  to  prune  the  existing  road  points  which  are  clearly  not  associated 
with  x*;. 

If  Xfc  belongs  to  an  exiting  road,  d\  follows  a  Chi-square  distribution  with  nx  degrees  of 
freedom,  where  nx  is  the  dimension  of  the  target  position  vector,  which  is  three  when  the 
road  is  represented  in  a  three-dimensional  space.  For  a  pre-specihed  probability  of  type  I 
error  a  ,  the  test  threshold  can  be  obtained  as 

Ta  =  F~2  (1  -  a)  (3) 

A.Tl'p 

where  FA,1  is  the  inverse  of  the  cumulative  distribution  function  of  a  Chi-square  distributed 
random  variable  with  nx  degrees  of  freedom.  When  d\  <  Ta,  we  can  decide  that  the  point 
belongs  to  the  existing  roads. 

2.2.2  Likelihood  with  Diffuse  Prior 

In  [11],  under  the  assumption  that  the  prior  of  the  true  target  state  x  is  a  Gaussian  distri¬ 
bution  with  an  arbitrary  mean  //  and  scaled  covariance  rr2P,  Kaplan  et  al.  showed  that  as 
cr2  — y  oo,  the  likelihood  of  two  tracks,  x*,  and  x/,  corresponding  to  the  same  target  state  x, 
becomes 

L(xfc,  X/)  = - - - -  exp  —  ^Ay  (Pfc  +  P i  —  P ki  ~  P ikV1  Akl  (4) 

|27r(Pfc  +  P,-PfcZ-PZfc)|2  L  2  J 

Therefore,  a  proper  test  statistic  is  obtained  in  a  straightforward  manner  as  follows 

d2  =  min  [log  |27r  (Pfc  +  P;  -  P ki  ~  Pa-)  I  +  (pk  +  P i~  P ki  ~  P/fc)_1  Afc^]  (5) 

Comparing  d2  in  (5)  with  d\  in  (2),  it  is  clear  that  they  are  almost  the  same,  except  that 
in  d2,  an  extra  penalty  term  has  been  introduced  to  favor  those  associations  resulting  smaller 
covariance  matrix  (Pfc  +  P/  —  Pfc;  —  P;*,). 
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2.2.3  Generalized  Likelihood 


If  the  target  state  estimate  can  be  treated  as  a  deterministic  vector,  one  can  use  the  gener¬ 
alized  likelihood  approach  [11],  Given  the  target  state  x  is  known,  the  likelihood  of  x*.  and 
Xi  is 


L(xfc,Xj|x)  =  - : 

|27rE|  ■ 

where  x2  =  [x£  xf]r,  I2  =  [I Uxxnx  I n, 


exp 


i  <X2 


I2x)r  E  1  (x2  -  I2x) 


1T  and 


(6) 


Pfc  P  kl 
P  Ik  P/ 


Since  x  is  unknown,  in  a  generalized  likelihood,  it  is  replaced  by  its  ML  estimate  based 
on  Xfc  and  x*,  which  is  the  fused  state  estimate  xj.  It  is  well  known  that  the  fused  state  xj 

is 

x/  =  P /  (7) 

where 

P  /  =  (I^E-Up-1  (8) 

Therefore,  the  GL  is  given  as  follows 


L(xfc,xi|x) 


1 

- -exp 

|27rE| 2 

1 

- rexp 

1 27rE  1 2 


1  T  -1 

--(x2-I2  X/)  E  (x2-I2x/) 
(X2E_1X2  -xJP^Xj) 


As  a  result,  the  test  statistic  based  on  the  GL  is 

d3  =  min  [log  |27tE|  +  (x^S_1x2  —  xJP^x/)] 


(9) 


(10) 


3  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

3.1  Conversion  between  Different  Coordinate  Systems  and  Dataset  Visualiza¬ 
tion 

Matlab  codes  are  developed  for  the  conversion  from  geographic  coordinates  to  topocentric 
coordinates,  and  for  the  display  of  the  road  map  information. 

In  order  to  process  the  synthetic  dataset,  it  is  required  to  transform  the  coordinates  of 
a  target  from  the  geographic  coordinate  system  to  the  topocentric  coordinate  system.  The 
earth  is  assumed  to  be  an  ellipsoid,  and  the  geographic  coordinates  consist  of  latitude  0, 
longitude  A,  and  the  ellipsoid  height  h  of  the  target.  A  topocentric  coordinate  system  is  a 
3-D  Cartesian  system  having  mutually  perpendicular  axes  U,  V,  W  with  an  origin  on  or 
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near  the  surface  of  the  Earth.  The  U-axis  is  locally  east,  the  V-axis  is  locally  north  and  the 
W-axis  is  up  forming  a  right-handed  coordinate  system. 

For  a  given  topocentric  origin  with  the  geographic  coordinates  0o,  A0,  and  ho,  the  conver¬ 
sion  of  a  target’s  geographic  coordinates  (</>,  A,  h )  to  its  topocentric  coordinates  (U,  V,  W), 
or  the  reverse  formula  to  convert  topocentric  coordinates  (U,  V. ,  W)  into  latitude,  longitude 
and  ellipsoidal  height  (</>,  A,  h),  has  been  discussed  in  details  in  [12], 

Based  on  the  formulas  described  in  [12] ,  Matlab  codes  for  the  conversion  from  geographic 
coordinates  to  topocentric  coordinates  have  been  developed.  Specifically,  the  earth  is  as¬ 
sumed  to  be  an  ellipsoid  with  the  equatorial  radius  of  6378137.0  meters,  and  the  polar 
radius  of  6356752.314  meters.  The  Salt  Flats  of  Utah  are  assumed  to  be  1286  meters  above 
the  sea  level.  The  topocentric  origin  is  taken  as  the  center  of  the  10km  x  10km  area,  with 
its  latitude  0o  =  40.7  degrees,  longitude  Ao  =  —113.875  degrees,  and  the  ellipsoid  height 
h0  =  1286  meters. 

The  truth  data  in  SyntheticData.mat  include  the  road  segments,  which  are  indicated  by 
the  geographic  coordinates  (latitude  and  longitude)  of  the  starting  and  end  way  points.  The 
Matlab  program  takes  the  truth  data,  converts  the  data  to  the  topocentric  coordinate  system, 
and  plots  the  road  network  in  the  topocentric  coordinate  system.  In  addition,  the  program 
also  converts  the  track  data  from  geographic  coordinate  system  to  topocentric  coordinate 
system.  In  Figure  1,  the  map  network  and  the  Erst  5  target  tracks  in  SyntheticData.mat 
have  been  plotted.  It  is  clear  that  the  tracks  are  very  close  to  the  corresponding  roads,  but 
they  do  not  coincide  with  each  other,  due  to  estimation  errors  in  the  track  data. 

3.2  Road  Conflation  Based  on  Radar  Positional  Data 

In  this  section,  an  efficient  mathematical  representation  of  the  road  is  proposed,  the  likelihood 
function  of  the  track  data  conditioned  on  the  road  parameters  is  derived,  and  solutions  to 
road  estimation  based  on  track  data,  track  to  road  association,  road  to  road  association,  and 
road  fusion  are  developed. 

3.2.1  Road  Representation 

In  the  proposed  approach,  roads  are  represented  by  piece-wise  linear  segments.  More  par¬ 
ticularly,  a  road  is  represented  by  the  following  straight  line  in  the  2-D  space: 

y  =  ax  +  b  i£[(s  fe]  (11) 

where  a  is  the  slope,  b  the  y-intercept.  and  £e  are  the  x-coordinates  of  the  start  and 
end  points  of  the  road,  respectively.  The  problem  of  road  estimation  is  to  estimate  the 
parameters,  including  a,  b,  £s,  and  £e,  based  on  the  radar  tracks. 
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Figure  1:  Road  network  and  the  first  five  tracks  in  SyntheticData.mat.  Solid  lines:  roads; 
circles:  tracks. 


3.2.2  Likelihood  of  Track  Data 

Let  us  assume  that  a  segment  of  track,  which  corresponds  to  the  same  road,  consists  of  N 
data  points: 


Z1  :N 


[x  1)  V\i  i  XN,  Vn\ 


J1 ) 


TiT 

lN1 


(12) 


where  z ;  =  [£*,  j/*]r  is  the  estimate  of  the  x  and  y  coordinates  of  the  target  at  time  i,  whose 
estimation  error  covariance  is  denoted  as  Xj.  Further,  the  parameters  of  the  estimation 
problem  are  denoted  as  s  =  [xi,  x2,  •  •  • ,  xjv,  a,  b]T ■  With  these  notations  and  assumptions, 
the  likelihood  of  z1:N  can  be  derived  as 


N 

p{  z1:jv|s)  =  Yl 
2=1 


1  (2i~Mi)TSi  1(zi-^i) 

- P  2 


where  /x?:  =  [xi  aXi  +  b}1 .  The  corresponding  log-likelihood  is 


N 

logp(zi;JV|s)  =  ^ 
2=1 


z  i 


I  R 


(13) 


(14) 


where  c  is  a  constant  independent  of  s. 
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3.2.3  Maximum  Likelihood  Estimation  of  Road 

The  maximum  likelihood  estimator  (MLE)  of  s  is  therefore 

N 

s  =  arg  min  VVzi  -  Aq)TE -  /**)  (15) 

S  L ' 

2=1 

We  derive  the  MLE  of  s  and  summarize  it  in  the  following  theorem. 

Theorem  1  The  optimal  MLE  of  s  based  on  position  information  in  radar  tracks  is  a  two- 
step  procedure.  First,  the  estimate  of  Xi,  denoted  as  fi,  should  satisfy  the  following  condition 


where  c  =  [1  a]T  and  dj 
rj  =  [a  b]T  is 

f] 


cTE“1di 


Zi  = 


V  i 


(16) 


cTE.-‘c 

[—Xi  b  —  y.i]T .  Plugging  (16)  back  into  (15),  the  MLE  of 


N 


arg  min  )  dj  E,  *dj  — 

71  <  ^ 


-Tsr'd,); 


2=1 


ctE7ic 


(17) 


Proof:  First,  we  re-write  the  cost  function  in  (15)  as 

N 


N 


2=1 


i  -  Zi)T^i  -  Zi)  =  ( CXi  +  di)T  1  ( CXi  +  d*) 

2=1 
N 

—  ^2  ‘'C  +  2crE,-1d ,Xi  +  d^Epd, 


(18) 


2=1 


where  the  first  equality  follows  from  the  identity 


A h  -  z i  =  [xi  -  Xi  axi  +  b-  yi]T 

=  [1  a}TXi  +  [-Xi  b  -  yi]T 

=  CXi  +  dj 


(19) 


Since  the  it\i  element  in  the  summation  in  (18)  is  only  a  quadratic  function  of  xV}  the 
optimal  estimator  of  x%,  must  satisfy  the  following  condition 


ti  = 


ctE71c 


(20) 


Substituting  this  condition  back  into  (18),  we  finally  have 

Wd, 


2=1 


cTE 


-l. 
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(21) 


Therefore,  the  MLE  of  77  is  (17),  and  the  MLE  of  ay  is 


ci.'nfZj  Mj(^) 

0(77)^0(77) 


(22) 

Q.E.D. 


Note  that  Theorem  1  simplifies  the  MLE,  an  optimization  problem,  significantly  when 
N  is  large.  The  original  MLE  problem  as  defined  in  (15)  involves  optimization  in  a  (iV  +  2)- 
dimensional  search  space,  whereas  the  optimization  problem  in  Theorem  1  requires  a  search 
in  merely  a  2-dimensional  space. 

We  also  derive  the  Cramer-Rao  lower  bound  (CRLB)  on  the  estimation  error  covariance 
matrix  for  any  unbiased  estimator  of  77,  which  is  provided  in  the  following  theorem. 


Theorem  2  The  Fisher  information  matrix  F,  which  is  the  inverse  of  the  CRLB  matrix, 
for  estimating  77  is 

N 

F  =  ^E,-G?.  (23) 

i=l 


where 


E,  = 


ir  i 


efE-'e 

efEr1  f  f'E-Jf 


el  srAf 

fTV-L 


G;  = 


cTE“1( 


(VE-‘e.)2  (c’-E.r-e.)  (c^pf) 

cTS-'e()  (c^E-'f)  (c’-E-'f)2 


=  [0  Xi]T ,  and  f  =  [0  1]J 


(24) 

(25) 


Proof:  Let  us  denote  L 
estimating  s  is 


Using  (14),  we  have 


logp(zi:7v|s),  then  the  Fisher  information  matrix  (FIM)  J  for 


J=E{V,£Vft| 


(26) 
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Let  us  define  Xi  =  ay  —  ay  and  iy  =  a  ay  +  b  —  fa.  Then,  it  is  easy  to  show  that  [ay  yi\T  has 
zero  mean  and  covariance  matrix  Ej,  and  (27)  can  be  re-written  as 

-[1  a]  S)"1  %  yif  =  -cT  E"1  [x{  UT  (28) 

Following  a  very  similar  procedure,  we  can  prove  the  following 


IT  N  N 

-  Y 1° s?rl  y^T  =  ~Y1  e?Xrl 


da 


■T'i  Vi 


i=  1 


i= 1 


and 


Now  we  have 


dL 

~db 


N 


N 


-  Y  [°  1]  Si  1  [Xi  Vi]1  =  -f  Y  Si  1  [®i  Vi. 


i=  1 


2—1 


E^lr)  >  =  E  E*r‘ [ii  iiYlii  ii]  E-“c 


Similarly,  we  can  prove  that 


^  ,  dL  dL  .  w  , 

E  w—  ^  =  0  V  1  ^  j 

dxi  dxj 


E 


dL  dL 
dxi  da 


=  cT  Ej 1ei 


(29) 

(30) 

(31) 

(32) 

(33) 


r  an 

\  day  db  J 


cT  Er'f 


(34) 


£ 


2—  1 


(35) 


and 


(36) 
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In  summary,  we  have 


J  =  £{VsLVfL} 


cT  Ej 

7c 

0 

0 

T 

C 

sr'ei 

T 

C 

E^f 

*  O 

cT  S2  7c  • 

0 

T 

C 

S2  1e2 

T 

C 

S2  7f 

0 

0 

•  cTE^c 

T 

C 

E  NxeN 

T 

C 

cT  Eb 

Xei 

cT  S2 1e2  • 

•  cT  E^ejv 

1^1=1 

ef  E"1^ 

2^i= 1 

ef  Er 

7f 

cT  Ef 

7f 

cT  S2-'f  • 

•  cT  E^f 

2^i=i 

ej  Er‘f 

2^i= 1 

fT  E- 

7f 

A  B 
Br  D 


where 


and 


A  = 


B 1 


cT  E,  7c 
0 

0 

14 


Ei  f 


c 

c 


Ej7c 


T  E2-1e2 


T  Eo  xf 


0 

0 


zT  E"1 


N 

N±eN 
N 


zT  E"1 


zT  Erdf 


D  = 


Vjv  pt 

The  CRLB  matrix  Cc  for  estimating  s  is 


Ef=1  ef  Er^ 

\N 


Ef=1ef  S“1f 

Lf 


Cc  =  J1 


A 


(38) 


(39) 


(40) 

(41) 


(42) 


Using  blockwise  matrix  inversion  [13],  we  know  that  the  lower  right  2x2  sub-matrix  cor¬ 
responding  to  the  estimation  of  77  =  [a  b]T  is  (D  —  B7  A  :B)  .  As  a  result,  the  FIM  for 

estimating  77  is  simply 

F  =  D  -  BTA  XB  (43) 

Plugging  (39),  (40),  and  (41)  into  (43),  we  can  finally  prove  the  theorem.  Q.E.D. 


From  Theorem  2,  it  is  clear  that  complexity  in  calculating  F  is  significantly  reduced  when 
N  is  large.  The  calculation  of  F  involves  the  inverse  of  matrices  with  smaller  dimensions 
(2x2),  instead  of  the  inverse  of  the  original  (N  +  2)  x  (IV  +  2)  FIM  J  for  s. 


3.2.4  Track  to  Road  Association 

To  test  whether  or  not  a  track  data  point  z  =  [x  y]T  with  covariance  matrix  E  belongs  to  a 
particular  road  with  parameter  estimate  77  =  [a  b]T  and  covariance  matrix  E^,  a  Chi-square 
test  is  adopted. 
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First,  denote  Ai  —  y  —  (ax  +  b ).  Then,  it  can  be  shown  that 

E[  A,]  =  0 

Var[Ai]  =  gTS  g  +  h%h  +  E,(l,  1)  x  E(l,  1) 


(44) 

(45) 


A2 

As  a  result,  the  test  statistic  is  set  as  yaj*A  - ,  which  given  a  Gaussian  assumption,  follows 

A2 

a  Chi-square  distribution  with  1  degree  of  freedom.  Hence,  if  yar*A  ^  <  H,  then  z  belongs  to 
the  road  parameterized  by  f),  where  t\  is  a  pre-defined  threshold  which  can  be  determined 
by  a  pre-specified  probability  of  type  I  error. 


3.2.5  Road  to  Road  Association 

To  test  whether  or  not  a  road  estimate  r)1  with  covariance  matrix  Zb  can  be  associated 
with  another  road  estimate  r)2  with  covariance  matrix  Z2,  a  similar  Chi-square  test  can  be 
adopted. 

First,  denote  A2  =  r]l  —  f]2-  Then,  it  can  be  shown  that 

E[  A2]  =  0  (46) 

Cov[A2]  =  Si  +  S,  (47) 

As  a  result,  the  test  statistic  is  set  as  A^(Ei  +  Z2)_1A2,  which  given  a  Gaussian  as¬ 
sumption,  follows  a  Chi-square  distribution  with  2  degrees  of  freedom.  Hence,  if  A^Zb  + 
Z2)_1A2  <  f2,  then  these  two  road  estimates  can  be  taken  as  to  correspond  to  the  same 
road,  t-2  is  a  pre-defined  threshold  which  can  be  determined  by  a  pre-specified  probability  of 
type  I  error. 

3.2.6  Road  Fusion 

If  a  road  estimate  r]l  with  covariance  matrix  Zi  can  be  associated  with  another  road  estimate 
f]2  with  covariance  matrix  Z2,  they  should  be  fused  to  get  a  more  accurate  road  estimate. 
The  fused  road  estimate  is 

f)  =  Z2(Zi  +  Z2)  lf]l  +  Zi(Zi  +  Z2)  1r)2  (48) 

and  the  corresponding  covariance  matrix  of  the  fused  estimate  is 

Z  =  Z2(Z1  +  Z2)-1Z1  (49) 

3.3  Road  Conflation  Based  on  Both  Radar  Positional  and  Velocity  Data 

In  Section  3.2,  only  position  estimates  in  the  radar  tracks  were  used  for  road  network  es¬ 
timation.  The  velocity  estimates,  which  also  contain  valuable  information  about  the  roads 
(especially  that  about  the  directions  of  the  roads),  were  not  employed.  In  this  section,  a 
road  estimation  approach  is  developed  to  incorporate  both  position  and  velocity  information 
in  the  radar  track  data.  The  maximum  likelihood  estimator  (MLE)  of  a  road  based  on  radar 
tracks  and  its  corresponding  Cramer-Rao  lower  bound  matrix  are  derived. 
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3.3.1  Likelihood  of  Track  Data 


As  in  Section  3.2,  roads  are  represented  by  piece-wise  linear  segments,  which  are  modeled  by 
straight  lines  in  the  2-D  space  with  slope  a  and  y-intercept  b.  It  is  reasonable  to  assume  that 
the  target  moves  towards  the  same  direction  as  that  defined  by  the  current  road  segment, 
except  when  it  makes  turns  at  road  intersections.  Hence,  we  have 


Vy  OjVx 


(50) 


where  vx  and  vy  are  the  target  velocities  along  the  x  and  y  axes,  respectively. 

Let  us  assume  that  a  segment  of  track,  which  corresponds  to  the  same  road,  consists  of 
N  data  points: 


Z1 :7V 


[-b  i ,  2/i, 


2/1 1 


‘^Ni  VNi  vxn  i  ^VNi 


1 T 


J1  ) 


T  ]T 

*N. 


(51) 


where  z%  =  [£*,  y*,  vXi,  vyi]T  is  the  state  estimate  of  the  target  at  time  i,  which  consists  of  the 
estimates  of  the  target  position  and  velocity  along  x  and  y  axes,  respectively.  The  estimation 
error  covariance  matrix  is  denoted  as  E*.  Further,  the  parameters  of  the  estimation  problem 
are  denoted  as  s  =  [x\,  vXl,  X2,  vX2,---,  xn,  vXn,  a,  b]T.  With  these  notations  and 
assumptions,  the  likelihood  of  z1:N  can  be  derived  as 


N 

p(z1:jv|s)  =  Yl 

i=l 


1 

- 6  2 

|2ttE  ^ 


where  //,-  =  [xi,  aXi  +  b ,  vXi,  avXi]T.  The  corresponding  log-likelihood  is 


N 

10gp(Zl;Jv|s)  =  Y2 

i— 1 


z  i 


+  c 


(52) 


(53) 


where  c  is  a  constant  independent  of  s. 


3.3.2  Maximum  Likelihood  Estimation  of  Road 

The  Maximum  Likelihood  Estimation  (MLE)  of  s  is  therefore 

N 

s  =  arg  min  yVzi  -  -  /xj  (54) 

S  L ' 

i=  1 

Analogous  to  Theorem  1,  we  derive  a  simplified  version  of  the  MLE  as  summarized  in 
the  following  theorem. 

Theorem  3  The  optimal  MLE  is  a  two-step  procedure.  For  1  <  i  <  N ,  let  us  denote 
Xj  =  [xi,  vXi]T .  The  estimate  of  x,,  denoted  as  xi7  should  stratify  the  following  condition, 

x,;  =  -  (CtE“1C)_1  C^r'm,  Mi 
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(55) 


where 


C 


1  0 
a  0 
0  1 
0  a 


and 

m,:  =  [-Xi,  b-jji,  -vXi,  -%]T. 
Plugging  (55)  back  into  (54),  the  MLE  of  rj  =  [a,  b]T  is 

N 


r]  =  arg  nun 
v 


£  [mJX-'m,  -  mfE-'C  (CTE,-1C)-1  CTE-‘m, 
2=1 


(56) 


The  proof  of  Theorem  3  is  similar  to  that  of  Theorem  1,  and  is  skipped  here  for  brevity. 
Again,  Theorem  3  simplifies  the  MLE  problem  significantly  when  N  is  large.  The  original 
MLE  problem  as  defined  in  (54)  involves  optimization  in  a  (2A^+2)-dimensional  search  space, 
whereas  the  optimization  problem  in  Theorem  3  requires  a  search  in  merely  a  2-dimensional 
space. 

The  Fisher  matrix  F  for  estimating  s  is 


F 

1  S 


where 


are  2x2  sub-matrices,  and 


1 

> 

0 

0 

1 - 

< 

•  o 

A2 

0 

A.2  ,n 

..  o 

••  0 

an 

p- 

..  ^ 
< 

L 

AT 

^2,77 

...  at 

L1N,r) 

- 1 

C* 

< 

Ai  =  CtE“1C 
A,;,,  =  CTE“1D.j 

N 

Xj  —  £  DfEpD, 
2—1 

0  0  " 

D,=  Xi  1 

0  0 

Vx,  0 


(57) 


(58) 

(59) 

(60) 

(61) 


The  CRLB  matrix  for  estimating  r]  —  [a,  6]r  can  be  obtained  by  first  taking  the  inverse 
of  the  (21V +  2)  x(2jV  +  2)  F,  matrix,  and  then  taking  the  lower-right  2x2  sub-matrix  of  F”1. 
However,  this  is  not  a  very  efficient  approach  especially  when  N  is  large  and  the  inversion  of 
a  large  Fs  is  involved.  As  we  did  in  the  proof  of  Theorem  2,  taking  advantage  of  the  special 
structure  of  Fs,  and  using  blockwise  matrix  inversion  [13],  a  much  more  efficient  approach  to 
calculate  the  Fisher  information  matrix  F,?  is  derived  and  provided  in  the  following  theorem. 
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Theorem  4  The  Fisher  information  matrix  for  estimating  q  based  on  both  position  and 
velocity  information  in  radar  tracks  is 

N 

F„  =  £  [DfS-‘D,  -  DfE-1C(CTE-'C)-1CTE-1Di]  (62) 

i= 1 

It  is  clear  that  the  calculation  of  F^  only  involves  manipulation  of  matrices  with  smaller 
dimensions  (such  as  2  x  2,  4  x  2,  2  x  4,  and  4  x  4.),  instead  of  the  original  (2IV  +  2)  x  (2A  +  2) 
matrix.  Now,  the  CRLB  matrix  for  estimating  ry  can  be  readily  obtained  by  taking  the 
inverse  of  F,;. 

3.3.3  Track  to  Road  Association  Using  both  Position  and  Velocity  Data 

To  associate  a  track  data  point  z  =  [x  y  vx  vy]T  with  covariance  matrix  £  to  a  road  with 
parameter  estimate  f]  =  [o  6] 1  and  covariance  matrix  a  Chi-square  test  is  adopted. 
Assuming  that  z  belongs  to  the  road,  we  denote  A2  =  [y  —  (ax  +  b)  vy  —  avx]T.  Then,  it  can 
be  shown  that 

E[  A2]  =  0  (63) 

where  0  is  a  all-zero  vector.  The  covariance  matrix  of  A2  is 

SA  =  E{A2A^} 

'  fT£f  +  gTE,g  +  £(l,l)E,?(l,l)  fTSh  +  gTE;/k  +  E(l, 3)Sr?(l,  1) 

_  fTEh  +  gTE})k  +  E(l,  3)E„(1, 1)  hTEh  +  kTS,;k  +  E(3,  3)E„(1, 1) 

in  which 

f  =  [—a  1  0  0]T,  h  =  [0  0  —  a  1]T,  g  =  [— x  —  1]T,  k  =  \—vx  0]T.  (65) 

The  test  statistic,  t  =  A^E^XA2,  follows  a  Chi-square  distribution  with  2  degrees  of 

freedom.  If  t  <  t2,  then  z  belongs  to  the  road  parameterized  by  f),  where  f2  is  a  pre-dehned 
threshold  corresponding  to  a  pre-specihed  probability  of  type  I  error. 

3.4  Overview  of  the  Road  Extraction  Algorithm 

Here  we  provide  an  overview  of  the  road  extract  algorithm.  Radar  tracks  are  used  as  building 
blocks  to  generate  road  estimates.  This  is  an  iterative  algorithm  where  a  single  radar  trajec¬ 
tory  is  processed  at  each  iteration.  During  each  iteration,  new  road  estimates  will  be  either 
fused  with  exiting  road  estimates  or  added  as  new  road.  As  more  tracks  are  integrated  into 
an  initially  empty  map,  the  accuracy  and  richness  of  the  map  improve.  A  post  processing 
algorithm  can  be  applied  to  merge  fragmented  roads  and  prune  isolated  false  roads.  The 
overview  of  the  algorithm  is  illustrated  in  Figure  2. 
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Figure  2:  Road  map  extraction  based  on  radar  tracks. 

4  RESULTS  AND  DISCUSSION 

4.1  Numerical  Examples  based  on  Position  Data  Only 

Based  on  the  theoretical  results  derived  in  Section  3.2,  we  develop  data  structure  for  road 
representation,  and  Matlab  codes  for  road  estimation  based  on  simulated  radar  tracks,  track 
to  road  association,  road  to  road  association,  and  road  fusion. 

In  Figure  3,  a  road  estimation/extraction  example  is  shown.  The  extracted  roads  are 
based  on  the  first  ten  radar  tracks.  As  shown  in  the  figure,  most  of  the  time,  the  estimated 
roads  are  reasonably  accurate  and  close  to  the  underlying  true  road  networks.  However,  for 
short  road  segments,  the  estimates  sometimes  are  not  very  accurate.  This  is  especially  clear 
for  part  for  the  polygon  circle  in  the  lower-left  corner  of  the  figure. 

In  Figure  4,  the  road  estimates  based  on  all  the  radar  tracks  are  presented.  As  shown  in 
the  figure,  the  estimated  roads  are  most  concentrated  around  the  true  roads.  But  there  are 
a  lot  of  fragmented  roads  and  some  false  roads.  To  alleviate  this  problem,  a  post  processing 
algorithm  is  applied  to  merge  fragmented  road  estimates  and  eliminate  isolated  false  road 
estimates.  Road  estimates  after  post  processing  (merging  and  pruning)  are  shown  in  Figure 
5.  It  is  clear  that  even  after  the  post  processing,  the  road  estimates  are  still  not  ideal,  and 
there  still  exist  fragmented  road  estimates  and  false  road  estimates. 

4.2  Numerical  Examples  Based  on  Both  Position  and  Velocity  Data 

In  Figure  6,  a  road  estimation/extraction  example  is  shown  based  on  the  first  ten  radar 
tracks,  where  both  position  and  velocity  information  is  used.  As  shown  in  the  figure,  most 
of  the  time,  the  estimated  roads  are  reasonably  accurate  and  close  to  the  underlying  trne 
road  networks.  However,  compared  to  Figure  3,  there  is  an  extra  false  road. 

In  Figure  7,  the  road  estimates  based  on  all  the  radar  tracks  are  presented.  As  shown 
in  the  figure,  the  estimated  roads  are  most  concentrated  around  the  trne  roads.  But  there 
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Figure  3:  Estimated  roads  based  on  the  first  ten  radar  tracks.  Blue  solid  lines:  true  roads; 
red  dash-dot  lines+circles:  estimated  roads. 
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Figure  4:  Estimated  roads  based  on  all  radar  tracks.  Blue  solid  lines:  true  roads;  red 
dash-dot  lines+circles:  estimated  roads. 
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Figure  5:  Estimated  roads  after  pruning  isolated  road  estimates.  Blue  solid  lines:  true  roads; 
red  dash-dot  lines-)- circles:  estimated  roads. 


Figure  6:  Estimated  roads  based  on  the  first  ten  radar  tracks.  Blue  solid  lines:  true  roads; 
red  dash-dot  lines+circles:  estimated  roads. 


Approved  for  Public  Release;  Distribution  Unlimited. 

17 


2000 

—  1000 

s 

0 

£ 

+J 

M  -1000 

o 

IS 

-2000 

-3000 


-5000 


Figure  7:  Estimated  roads  based  on  all  radar  tracks.  Blue  solid  lines:  true  roads;  red 
dash-dot  lines+circles:  estimated  roads. 

are  a  lot  of  fragmented  roads  and  some  false  roads.  To  alleviate  this  problem,  again  a  post 
processing  algorithm  is  applied  to  merge  fragmented  road  estimates  and  eliminate  isolated 
false  road  estimates.  Road  estimates  after  post  processing  (merging  and  pruning)  are  shown 
in  Figure  8.  It  is  clear  that  even  after  the  post  processing,  the  road  estimates  are  still  not 
ideal,  and  there  still  exist  fragmented  road  estimates  and  false  road  estimates.  However, 
comparing  Figure  5  and  Figure  8,  the  estimation  performance  improvement  achieved  by 
using  both  position  and  velocity  information  is  very  clear. 

4.3  Initialization  using  velocity  data 

From  the  experiments,  it  is  found  that  the  initialization  of  road  estimate  is  very  important 
and  has  direct  impacts  on  the  road  estimation  performance.  So  far  in  the  initialization  pro¬ 
cess,  only  radar  track’s  position  information  has  been  used  to  estimate  the  road  parameters 
a  and  b.  Since  the  radar  track’s  velocity  data  contain  the  road  direction  information,  in¬ 
corporation  of  this  information  has  the  potential  to  improve  the  road  initialization,  and  to 
improve  the  following  road  estimation  process.  The  following  modification  has  been  made 
to  include  the  velocity  information.  Let  us  define  the  target’s  heading  at  time  i  as 

Oii  =  atan2 (vy.,  vXi) 

where  vXi  and  %  are  the  target’s  velocity  estimates  along  x  and  y  axes  at  time  i  respectively. 
Two  consecutive  radar  track  data  points  at  time  i  and  i  +  1  are  used  to  initialize  a  road 
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Figure  8:  Estimated  roads  after  pruning  isolated  road  estimates.  Blue  solid  lines:  true  roads; 
red  dash-dot  lines+circles:  estimated  roads. 

estimate  only  if  the  two  directions  CKj  and  Qj+1  are  very  close  to  each  other.  Mathematically, 
the  condition  is  |a,  —  cp+i|  <  77,  or  (cp  +  7r  —  al+\  \  <  77,  or  |a;  —  7r  —  a!j+i|  <  77,  where  77  is  a 
pre-specified  threshold. 

In  this  experiment,  we  set  77  as  5°.  In  Figures  9  and  10,  the  road  estimates  which  are 
initialized  without  velocity  information  and  those  initialized  with  the  velocity  information 
are  compared.  It  is  clear  that  the  latter  provides  us  better  road  estimation  solutions  with 
less  false  roads. 

4.4  Road  Consistence  Test 

To  further  reduce  the  false  road  estimates,  a  road  estimate  consistence  test  is  proposed.  Once 
a  new  road  estimate  is  generated,  all  the  track  data,  which  have  been  used  to  generate  the 
road  estimate,  are  subject  to  a  track  to  road  association  procedure,  as  discussed  in  Section 
3.3.3.  If  all  the  track  data  points  passed  the  test,  meaning  that  they  are  associated  to  the 
new  road  estimate,  then  the  new  road  estimate  is  accepted  and  used  to  update  the  global 
road  estimates;  otherwise,  the  new  road  estimate  does  not  have  a  good  “match”  with  the 
track  data  and  it  is  discarded.  In  this  way,  the  possible  false  road  estimates  are  removed 
from  the  road  database. 

In  Figures  11  and  12,  the  algorithm  involving  road  consistence  tests  is  compared  against 
that  without  the  consistence  test.  It  is  clear  that  the  former  has  much  better  performance 
with  more  consistent  road  estimates. 
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Figure  9:  Estimated  roads  using  the  first  20  radar  tracks.  Roads  are  initialized  without 
velocity  data.  Blue  solid  lines:  true  roads;  red  dash-dot  lincs+circles:  estimated  roads. 


Figure  10:  Estimated  roads  using  the  first  20  radar  tracks.  Roads  are  initialized  with  velocity 
data.  Blue  solid  lines:  true  roads;  red  dash-dot  lines+circles:  estimated  roads. 


Approved  for  Public  Release;  Distribution  Unlimited. 
20 


4000 


t 

o 


J3 


— ©-© 

-4000  -3000  -2000  -1000  0  1000  2000  3000  4000  5000 

East  (m) 


Figure  11:  Estimated  roads  using  the  first  100  radar  tracks  without  road  consistence  tests. 
A  roads  are  initialized  without  velocity  data.  Blue  solid  lines:  true  roads;  red  dash-dot 
lines+circles:  estimated  roads. 


Figure  12:  Estimated  roads  using  the  Erst  100  radar  tracks  with  road  consistence  tests.  Blue 
solid  lines:  true  roads;  red  dash-dot  lines+circles:  estimated  roads. 
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4.5  Road  Merging 

As  we  can  see  in  Fig.  12,  the  road  estimates  are  fragmented,  with  many  small  road  segments, 
some  of  which  really  belong  to  single  roads.  To  combine  the  fragmented  road  estimates  that 
belong  to  the  same  road,  Matlab  codes  have  been  developed  for  a  road  merging  algorithm. 
The  algorithm  works  as  follows. 

The  input  of  the  algorithm  is  the  road  estimates  based  on  the  radar  data,  and  the  output 
is  the  merged  road  estimates.  Two  sets  of  road  estimates  are  maintained.  The  first  set 
contains  the  merged  road  estimates,  which  is  empty  initially.  The  second  set  includes  the 
unprocessed  road  estimates.  At  each  iteration,  one  unprocessed  road  estimate  from  the 
second  set  is  compared  with  the  road  estimates  in  the  first  set.  If  it  can  be  associated 
with  a  merged  road  estimate  from  the  first  set,  using  a  threshold  t2\  as  discussed  in  Section 
3.2.5,  then  the  unprocessed  road  estimate  will  be  fused  with  the  corresponding  merged  road 
estimate,  as  described  in  Section  3.2.6,  and  it  will  be  deleted  from  the  second  set.  If  it  can 
not  be  associated  with  any  road  in  the  first  set,  then  road  to  road  association  tests  will 
be  performed  to  see  if  it  can  be  associated  with  any  road  in  the  second  set  using  another 
threshold  f22  •  If  the  unprocessed  road  is  associated  to  another  unprocessed  roads,  these  two 
roads  will  be  fused,  the  fused  road  will  be  added  to  the  first  set,  and  both  of  the  original 
roads  will  be  deleted  from  the  second  set.  If  the  unprocessed  road  can  not  be  associated 
with  any  road  in  both  sets,  then  it  will  be  added  in  the  first  set  as  a  new  merged  road  and 
deleted  from  the  second  set  as  an  unprocessed  road. 

Based  on  the  first  1000  radar  tracks,  258  road  estimates  have  been  generated  and  plotted 
in  Fig.  13.  The  road  merging  algorithm  has  been  applied  to  the  258  road  estimates,  and  the 
output,  the  merged  roads,  really  depends  on  the  thresholds  t2i  and  t22  used  in  the  algorithm. 
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Figure  13:  Estimated  roads  using  the  first  1000  radar  tracks  with  road  consistence  tests. 
Blue  solid  lines:  true  roads;  red  dash-dot  lines+circles:  estimated  roads. 
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Ill  Table  1,  the  number  of  merged  road  estimates  is  listed  as  a  function  of  the  thresholds. 
Based  on  the  numerical  results,  it  is  clear  that  the  threshold  h\  and  h 2  play  a  crucial  role 
in  road  merging  performance.  A  larger  threshold  means  that  it  is  easier  to  associate  a  road 
to  another  road,  and  hence  it  results  in  a  smaller  number  of  merged  road  estimates. 


Table  1:  Number  of  merged  roads  (Nm)  based  on  258  unprocessed  road  estimates 


hi  —  h2 

13.82 

23.03 

32.24 

41.45 

50.66 

59.87 

69.14 

Nm 

206 

173 

142 

132 

121 

116 

110 

When  t-2i  =  t‘22  =  69.14  are  used,  a  total  of  110  merged  roads  have  been  obtained  and 
plotted  in  Fig.  14.  Compared  to  Fig.  13,  in  this  figure,  there  is  less  number  of  road  estimates 
due  to  the  road  merging  procedure. 


Figure  14:  Merged  road  estimates.  Blue  solid  lines:  true  roads;  red  dash-dot  lines+circles: 
estimated  roads. 

4.6  Evaluation  of  Road  Estimation  Accuracy 

In  this  subsection,  we  evaluate  the  road  estimation  accuracy  of  the  proposed  algorithm  with 
different  amount  of  radar  track  data.  In  this  example,  the  22nd  road  as  shown  in  Figure  15  is 
chosen  to  illustrate  the  estimation  accuracy.  We  define  the  absolute  estimation  error  in  road 
parameters  a  and  b  as  \a  —  a\  and  \b  —  b\,  respectively,  and  compare  them  with  the  theoretical 


Approved  for  Public  Release;  Distribution  Unlimited. 

23 


Figure  15:  The  22nd  road,  denoted  by  diamonds. 


CRLB  on  the  standard  deviations  (s.d.s)  of  the  estimation  errors  for  a  and  b.  The  results 
based  on  different  number  of  radar  tracks  are  provided  in  Table  2.  From  this  table,  it  is 
clear  that  as  more  radar  track  data  are  fused  to  estimate  the  road,  the  road  estimation  error 
decreases,  and  estimation  accuracy  improves  significantly.  Further,  the  calculated  CRLB  is 
comparable  to  the  actual  estimation  error,  and  it  indeed  provides  a  very  useful  estimation 
performance  measure. 


Table  2:  Road  estimation  errors 


Number  of  radar  tracks 

10 

20 

50 

100 

|  a  —  a 

0.028 

0.0011 

0.0034 

0.0026 

CRLB  on  s.d.  in  a 

0.01 

0.0041 

0.003 

0.003 

\b~'b\ 

48.91 

6.86 

3.14 

3.04 

CRLB  on  s.d.  in  b 

18.7 

7.06 

5.87 

5.86 

5  CONCLUSION 

In  this  project,  we  proposed  a  framework  to  generate  accurate  road  network  map  based 
on  radar  tracks  of  ground  targets.  An  ML  road  estimator  has  been  developed  to  estimate 
roads  using  radar  tracks,  and  the  corresponding  CRLB  matrix  was  derived  as  well.  Track  to 
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road  correlation/association  approaches  have  been  proposed  to  associate  track  segments  to 
existing  road  segments.  To  enhance  the  association  performance,  statistics  and  information 
extracted  from  the  track  estimate,  including  both  position  and  velocity  estimates,  have  been 
used.  For  tracks  overlapping  with  any  existing  roads,  a  track  to  road  fusion  algorithm  was 
proposed  to  fuse  the  track  with  the  associated  road  segments  to  improve  the  road  map 
accuracy.  Some  numerical  examples  have  been  provided  to  demonstrate  the  effectiveness  of 
the  proposed  approaches. 
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List  of  Acronyms 


CRLB  Cramer-Rao  Lower  Bound 

DTW  Dynamical  Time  Warping 

GMTI  Ground  Moving  Target  Indicator 

ML  Maximum  Likelihood 

MLE  Maximum  Likelihood  Estimator 

NIM  A  National  Imagery  and  Mapping  Agency 

SAR  Synthetic  Aperture  Radar 

STAP  Space-Time  Adaptive  Processing 

USGS  United  States  Geological  Survey 
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